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Abstract. 

Recent progress in numerical relativity has enabled us to model the non- 
perturbative merger phase of the binary black-hole coalescence problem. Based on 
these results, we propose a phenomenological family of waveforms which can model 
the inspiral, merger, and ring-down stages of black hole coalescence. We also construct 
a template bank using this family of waveforms and discuss its implementation in the 
search for signatures of gravitational waves produced by black-hole coalescences in the 
data of ground-based interferometers. This template bank might enable us to extend 
the present inspiral searches to higher-mass binary black- hole systems, i.e., systems 
with total mass greater than about 80 solar masses, thereby increasing the reach of 
the current generation of ground-based detectors. 



1. Introduction 

The first generation of ground-based gravitational wave detectors [TJ [5J [3] are currently 
operating at unprecedented levels of sensitivity and the LIGO detectors, in particular, 
have attained their design goals over a broad frequency range. The data from these 
detectors has been used to search for a wide variety of gravitational-wave sources including 
coalescing binary black- hole systems (see e.g., [HE]). In parallel with these experimental 
and observational achievements, a series of breakthroughs has occurred in numerical 
simulations of binary black hole systems [HI O [5] . Long-term evolutions of inspiralling 
black holes that last for several orbits have been obtained with several independent 
codes [U [TOl HU H21 [El El [T5J [16] , and accurate gravitational- wave signals have been 
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computed. It is now possible, in principle, to use these numerical-relativity results in 
astrophysical searches for gravitational waves. However, the high computational cost of 
these simulations makes it unfeasible to numerically generate all the necessary waveforms 
to cover the parameter space that needs to be searched. It is therefore necessary, 
at the present time, to use results from post-Newtonian (PN) theory to extend the 
waveforms obtained from numerical relativity (NR). The issue of matching NR and 
PN waveforms for equal-mass binary black-hole systems, and using numerical relativity 
results in gravitational- wave searches has been considered previously [T7l [TBI HU [20] . We 
generalize this to unequal mass systems and suggest a phenomenological template bank 
parametrized only by the masses of the two individual black holes. This template bank 
could be used to search for binary black-hole signals in data from current and future 
generations of gravitational- wave detectors. Our phenomenological approach is motivated 
by the work of Buonanno et al |21j . but using only physical parameters and generalized 
to include recent results from numerical relativity. 

In this paper we combine restricted 3.5PN waveforms [22] with results from NR 
simulations to construct "hybrid" waveforms for the quasi-circular inspiral of non- 
spinning binaries with possibly unequal masses. Restricting ourselves to the leading- 
order quadrupole modes, we find that the hybrid waveforms can be approximated by 
phenomenological analytical waveforms with fitting factors > 0.99 in the total mass range 
between 30 and 130 Mq for Initial LIGO [23] . For our analysis we use numerical waveforms 
obtained with two independent codes: (i) long waveforms (12 cycles) from equal mass 
binaries, provided by the AEI-CCT groups, using their CCATIE code [24] based on the 
Cactus framework [25] and Carpet mesh-rcfincmcnt driver [26]; (ii) waveforms from an 
unequal mass parameter study presented by the Jena group in [27], which have been 
obtained with the BAM code [13j . An analysis of these waveforms focusing on ring-down 
and higher modes has been presented in [55] . 

The rest of this paper is organized as follows: Section [2] summarizes the numerical 
simulations and how numerical waveforms have been computed. In Section [3] hybrid 
waveforms are produced by matching PN and NR waveforms. In Section [4] we 
propose a family of phenomenological waveforms in the Fourier domain, and study their 
impact for detection and parameter estimation by computing the fitting factors of the 
phenomenological waveforms with the hybrid ones. We also parametrize the best-matched 
phenomenological waveforms in terms of the physical parameters. Section [5] shows the 
astrophysical range of a search using the full coalescence waveforms and a preliminary 
comparison with other searches. Finally, Section [6] concludes with a summary of our 
results and plans for future work. 

2. Numerical simulations 

Both the BAM [13] and CCATIE [24j codes are finite-difference mesh-refinement codes 
solving the Einstein equations within the "moving puncture" framework El ESI [16] . 

In the wave-zone, sufficiently far away from the source, the spacetime metric can 
be accurately described as a perturbation of a flat background metric; let h a b denote the 
metric perturbation where a, b denote spacetime indices. Let t be the time coordinate used 
in the numerical simulation to foliate the spacetime by spatial slices. Working, as usual, 
in the transverse-traceless (TT) gauge, all the information about the radiative degrees of 
freedom is contained in the spatial part hij of h a b, where i, j denote spatial indices. Let 
us use a coordinate system (x,y, z) on a spatial slice so that the z-axis is along the total 
angular momentum of the binary system at the starting time. Let i be the inclination 
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angle from the z-axis, and let 4> be the phase angle and r the radial distance coordinates 
so that (r, I, <fi) are standard spherical coordinates in the wave-zone. 

Working in the TT gauge, the radiative degrees of freedom in h a t> can be written, as 
usual, in terms of two polarizations h + and h x : 

hij = h + (e + )ij + h x (e x )ij, (1) 
where e +iX are the usual basis tensors for TT tensors in the wave frame 

= ktj - 4>i(j)j , and = U<t>j + tj4>i ■ (2) 

Here i and <p ar e the unit vectors in the i and 4> directions respectively. The wave, of 
course, propagates in the radial direction. 

In our numerical simulations, the gravitational waves arc extracted by two distinct 
methods. The first one uses the Newman-Penrose Weyl tensor component ^4 (see e.g. 
|30| ) which, in an appropriate gauge is a measure of the outgoing transverse gravitational 
radiation in an asymptotically flat spacetime. By measuring that the peeling property 
(whereby '5 4 falls off as 1/r) is satisfied, we have determined that the gauge we are using 
does provide a good approximation to within the accuracy required for this study. In the 
wavezone it can be written in terms of the complex strain h = h + — ih x as |31j . 

h = lim / dt' [ dt" * 4 - (3) 

r ^°°Jo Jo 

An alternative method for wave extraction, which has a long history in numerical relativity, 
determines the waveform via gauge-invariant perturbations of a background Schwarzschild 
spacetime, via the Zerilli-Moncrief formalism (see [32] for a review). In terms of the even 

(Qim ) an d 0( id {Qim ) Parity master functions, the gravitational wave strain amplitude is 
then given by 

^ E (QL ~ i f QL(t') dt') +o(±). (4) 



Results from the BAM code have used the Weyl tensor component ^4 and Eq. ([3]), 
with the implementation described in [13]. While the CCATIE code computes waveforms 
adopting both methods, the AEI-CCT waveforms used here were computed using the 
perturbative extraction and Eq. |4]). Beyond an appropriate extraction radius, the two 
methods for determining h are found to agree very well for moving-puncture black-hole 
evolutions of the type considered here [TB] ■ 

It is useful to discuss gravitational radiation fields in terms of spin- weighted ,s = — 2 
spherical harmonics Yf m , and in this paper we will only consider the dominant £ = 2, m = 
±2 modes (see [55] for the higher I contribution in the unequal mass case), with basis 
functions 

r2 - 22 "/iE (1 - cost)V2i "' y ^^( 1+cost ) 2e2i '- ( 5 ) 

Our "input" NR waveforms correspond to the projections 



h em = (F £m 2 ,h) = / d<t> \\Y~2 sintd/,, (6) 
Jo Jo 

of the complex strain h, where the bar denotes complex conjugation. In the cases 
considered here, we have equatorial symmetry so that U22 = ^2-2, and 

h(t) = \f^^ ((1 + cost) 2 Mf) + (1 - cost) 2 h 22 (i)) . (7) 

In practice, we choose l = 0, thus h(i) = 4w (122^) ~ 0.631 b22(t). 
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Figure 1. NR waveform (black) from an equal- mass simulation, along with the 'best- 
matched' 3.5PN waveform (red). The Hybrid waveform constructed from the above is 
also shown (dashed line). 



3. Matching post-Newtonian and numerical relativity waveforms 



Once the PN and NR waveforms are generated, we produce a set of hybrid waveforms by 
matching them in an overlapping time interval t\ < t < t 2 . The obvious assumption in 
this procedure is that such an overlapping region exists and that in it both approaches 
yield the correct waveforms. 

Each time-domain waveform h(t, fi) is parametrized by a vector \x = {M, r\, <po, to}, 
where M = mi + is the total-mass of the binary, 77 = m,\mil (jn\ +7712) 2 is the 
symmetric mass-ratio, is the initial phase and t is the start time of the waveform. We 
match the PN waveforms /i+ N x (t, /1) and the NR waveforms h^ R x (t, u) Q by minimizing 
the square difference between the respective polarizations, i.e., 



5 = 



:=+,* 



*2 



K [t, ti)-a \ (t, u) 



(1/ 



(8) 



The minimization is carried over the parameters \x of the PN waveform and an amplitude 
scaling factor a [§|. The hybrid waveforms are produced by combining the 'best-matched' 
PN waveforms and the NR waveforms in the following way: 



h +tX (t, no) 



ift<h 



a Q Th + M,u) + (l- T)h + x (t,Ha) if *i < * < *2 ( 9 ) 



NR , . 

a h + ^(t,u) 



where Ho and ciq denote the values of \x and a for which 5 is minimum, and r = 
(i — t\)/{t2 — t\) is a linearly- increasing weighting function, such that < r < 1. 

An example set of hybrid waveforms is shown in Fig. [TJ The numerical waveform 
(black line) from an equal-mass (77 = 0.25) simulation by the AEI-CCT group is matched 
with a 3.5PN inspiral waveform (red line) over the matching region 100A/ < t < 850M. 
The hybrid waveform (dashed line) is constructed by combining the above as per Eq.([9]). 

I The parameters v are taken from the same set as fx, but with different values. 

§ The amplitude scaling factor was introduced in order to accommodate possible errors in the amplitude 
of the NR waveforms due to, for example, the finite radius of the extraction sphere. The best-matched 
value of a, however, was found to be 1 ± 0.08. 
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The robustness of the matching procedure can be tested by computing the overlaps 
between hybrid waveforms constructed with different matching regions. If the overlaps are 
very high, this can be taken as an indication of the robustness of the matching procedure. 
A more detailed discussion of this will be presented in [33] . 

Fig. [2] shows the hybrid waveforms of different mass-ratios in the Fourier domain. 
In particular, the panel on the left shows the amplitude of the waveforms in the Fourier 
domain, while the panel on the right shows the phase. These waveforms are constructed 
by matching 3.5PN waveforms with the NR waveforms from the unequal mass (0.16 < 
rj < 0.25) simulations by the Jena group. In the next section, we try to parametrize these 
Fourier domain waveforms in terms of a set of phenomenological parameters. 

4. The phenomenological template bank 

An obvious issue when using the above constructed hybrid waveforms directly as detection 
templates is that it might be computationally very expensive to compute enough NR 
waveforms to cover the entire parameter space densely enough. In this section, we propose 
a phenomenological waveform family which has more than 99 % overlaps with the hybrid 
waveforms in the detection band of the Initial LIGO detectors. We also show how this 
phenomenological waveform family can be mapped to the physical parameters (M and n), 
so that the template bank, at the end, is two-dimensional. 

4-1. Phenomenological waveforms 

We write our phenomenological waveform in the Fourier domain as 

uW^AesWf?*'*™. (10) 




Figure 3. Fitting factors of the hybrid waveforms with the phenomenological 
waveform family. Horizontal axis shows the symmetric mass ratio of the binary, while 
different colours/markers correspond to different total masses. 



where A c ff(f) is the amplitude of the waveform in the frequency domain, which we choose 
to write in terms of a set of "amplitude parameters" a — {/mcrg, /ring, /cut, &} as 

(///mcrg) ^ if / < /mcrg 

Aff (/) = I (f/f mcrg y 2/3 if / mcrg < / < / ring (ii) 

> W £(/, /ring, cr) if / ring < / < /cut 

In the above expression, 

CiL AinS ' ^ ^ (/-/ringV^ ' (12) 

represents a Lorcntzian function of width a centered around /ring- The normalisation 
constant w is chosen in such a way that -A G ff(/) is continuous across the "transition" 
frequency / ri n g , ie., 

-2/3 

(13) 



^ / / 



/ ring 



2 \ /me 

where we use / cu t as the cutoff frequency of the template and /merg as the frequency at 
which the power-law changes from /~ 7 / 6 to f~ 2 ^ 3 (as noted previously in [17] for the 
equal-mass case). 

Taking our motivation from the stationary-phase expansion of the gravitational-wave 
phase, we write the effective phase ^ c s(f) as an expansion in powers of /. 

* cff (/) = 2tt /t + + i'o r 5/3 + i'2 r 1 + ^3 r 2/3 + ^ r 1/3 + v 6 / 1/3 , (14) 

where to is the time of arrival, 4>o 1S the frequency-domain phase offset, and (3 = 
{"00, ip2> "03, "04, "06} are the "phase parameters", that is the set of phenomenological 
parameters describing the phase of the waveform. 
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The fitting factors [34j of the hybrid waveforms with the family of phenomenological 
waveforms are shown in Fig. [3] over the parameter range of 30 < M/M@ < 130 and 
0.16 < rj < 0.25, using the Initial LIGO noise spectrum. It is quite apparent that the 
fitting factors are almost always greater than 0.99, thus underlining the effectiveness of 
the phenomenological waveforms in reproducing the hybrid ones. As an example, in Fig.[¥J 
we plot the hybrid waveform and the best-matched phenomenological waveform from the 
M = 40 Mq, r) = 0.25 binary. 
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Figure 4. Hybrid waveform h(t) and the best-matched phenomenological waveform 
u(t) in the time domain for a M = 40 Mq, fj = 0.25 binary system. u(t) is computed 
by taking the inverse Fourier transform of the phenomenological waveform «(/). Both 
waveforms are normalised with respect to the Initial LIGO noise spectrum. 



4-2. From phenomenological to physical parameters 

It is possible to reparametrise the best-matched phenomenological waveforms in terms 
of the physical parameters of the hybrid waveforms. In Fig. we plot the amplitude 
parameters c* m ax of the best-matched phenomenological waveforms against the physical 
parameters of the binary. Similarly, the phase parameters /3 max of the best-matched 
phenomenological waveforms are plotted against the physical parameters of the binary 
in Fig. [6] The linear polynomial fits to the data points serve as a numerical "look-up 
table" to go from the physical parameters {M, rj} to the best-matched phenomenological 
parameters {a m ax, Anax}- Thus the template bank lives on a two-dimensional manifold 
(parametrised by M and 77) embedded in a higher-dimensional space. 

It might be worth stressing that the search will be carried out over M and 77, and 
not over the phenomenological parameters. The phenomenological parameters ct and /3, 
are constrained by the numerical look-up tables (see Figs. O and [S]), and only serve as an 
intermediate step in generating the template waveforms. 

5. The astrophysical range and comparison with other searches 

The template family proposed in this paper can be used for coherently searching for all 
the three stages (inspiral, merger, and ring-down) of the binary coalescence, thus making 
this potentially more sensitive than searches which look at the three stages separately. 
Fig. [7] compares the sensitivity of the searches using different template families. What 
is plotted here is the distance at which an optimally-oriented, equal-mass binary would 




Figure 5. Best-matched amplitude parameters ct max in terms of the physical 
parameters of the binary. The horizontal axis shows the symmetric mass-ratio of 
the binary and different colors/markers correspond to different total masses. Linear 
polynomial fits to the data points are also shown. 




Figure 6. The same as in Fig. [5] but for the phase parameters (3, 
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produce an optimal signal-to-noise ratio (SNR) of 8 at the Initial LIGO noise spectrum. 
The dotted line corresponds to a search using PN templates truncated at the Schwarzschild 
innermost stable circular orbit (ISCO), the dashed line corresponds to a search [35] using 
ring-down templates, and the solid line to a search using all three stages of the binary 
coalescence using the template bank proposed here. The horizontal axis reports the total- 
mass of the binary, while the vertical axis the distance in Mpc. It is quite evident that, 
for binaries with 50 < M/Mq < 140, the "coherent search" using the new template family 
is considerably more sensitive than any other search considered here. In a forthcoming 
paper [33] , we will quantify and extend this comparison to include other approaches which 
model all the three stages of black hole coalescence such as, for example the Effective One 
Body approach [36] , or alternatively, combining the results from separate inspiral and ring 
down phases using a coincidence analysis (see e.g. [37] )■ It is worth pointing out that 
regardless of the sensitivity, a coherent template bank such as the one proposed here is 
likely to be technically easier to implement than a inspiral-ring down coincidence analysis. 
While the method proposed here is just another template bank search, the inspiral-ring 
down coincidence would require tuning the individual inspiral and ring down searches, as 
well as finding the appropriate coincidence windows and other pipeline parameters. 

However, while this looks promising, we emphasize that it is important to treat Fig.[7j 
as only a preliminary assessment; fitting factors are not the only consideration for a 
practical search strategy. It is also very important to consider issues which arise when 
dealing with real data. For example, false alarms produced by noise artifacts, might well 
determine the true sensitivity of the search, and these artifacts will inevitably be present in 
real data. This is however beyond the scope of the present work, and further investigation 
is required before we can properly assess the efficacy of our phenomenological template 
bank in real- life searches. 

6. Summary and outlook 

Making use of the recent results from numerical relativity we have proposed a 
phenomenological waveform family which can model the inspiral, merger and ring-down 
stages of binary black-hole coalescence. We first constructed a set of hybrid waveforms 
by matching the NR waveforms with the analytical PN waveforms. Then, we analytically 
constructed phenomenological waveforms which approximated the hybrid waveforms. The 
family of phenomenological waveforms that we propose was found to have fitting factors 
larger than 0.99 with the hybrid waveforms in the detection band of Initial LIGO. We have 
also shown how this phenomenological waveform family can be parametrized in terms of 
the physical parameters (M and 77) of the binary, so that the template bank is, at the end, 
two dimensional. This phenomenological waveform family can be used to densely cover 
the parameter space, avoiding the computational cost of generating numerical waveforms 
at every grid point in the parameter space. We have also compared the sensitivity of a 
search using this template family with other searches. This search might enable us to 
extend the mass-range of the present inspiral searches to higher mass (> 80 M Q ) systems. 
In the mass-range 50 Mq to 140 Mq, this search could be significantly more sensitive than 
the search using the standard PN inspiral templates and quasi-normal mode ring-down 
templates. 

The numerical "look-up tables" to go from the physical parameters to the 
phenomenological parameters (see Figs. [5] and [5]) can be replaced by analytical functions 
of M and r\. This makes it easier to compute the parameter space metric used for 
template placement [38] and will be studied in future work. Our plans for future work 
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Figure 7. A preliminary assessment of the phenomenological template bank. This 
figure plots, as a function of the total mass, the distance to an optimally-oriented, 
equal-mass binary which can produce an optimal SNR of 8 at the Initial LIGO noise 
spectrum. The dotted line corresponds to a search using PN templates truncated at 
ISCO, the dashed line corresponds to a search using ring-down templates, and the 
solid line to a search using the template family proposed in this paper. The ring 
down horizon distance is computed assuming that e = 0.7% of the black hole mass is 
radiated in the ring down stage. Since the value of this parameter has some amount 
of uncertainty in it, we have also included the shaded region in the plot corresponding 
to 0.18% < e < 2.7%. 



also include the study of the robustness of the matching procedure used to construct the 
hybrid waveforms by considering different matching regions and PN waveforms of different 
order, and of course, to eventually construct a realistic search pipeline which incorporates 
numerical relativity waveforms in gravitational wave searches. 
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